Statistical learning: support vector machines

Objectives

• Define the maximal margin classifier
• Define the support vector classifier and discuss the logic of this approach
• Define support vector machines (SVM) and non-linear decision boundaries
• Apply SVM classification to example data sets and compare with alternative statistical learning models

library(tidyverse)
library(forcats)
library(broom)
library(modelr)
library(tree)
library(randomForest)
library(stringr)
library(ISLR)
library(titanic)
library(rcfss)
library(pROC)
library(gbm)
library(e1071)
library(grid)
library(gridExtra)

options(digits = 3)
set.seed(1234)
theme_set(theme_minimal())

Support vector machines (SVMs) are a popular statistical learning method for classification tasks.¹ SVMs build on several important concepts, that while related are distinct from one another. We will first discuss the logic of these individual components, then demonstrate how to estimate and interpret SVMs, and compare model results using this method to other statistical learning procedures we have discussed so far.

Maximal margin classifier

Hyperplanes

In \(p\)-dimensional space, a hyperplane is a flat subspace of \(p - 1\) dimensions that is affine (does not need to pass through the origin). In two dimensions, a hyperplane is a flat one-dimensional subspace (also known as a line). In three dimensions, a hyper plane is a flat two-dimensional subspace (also known as a plane). In higher dimensions it gets harder to visualize this concept, but the definition still holds true.

In two dimensions, the mathematical equation for a hyperplane is:

\[\beta_0 + \beta_1 X_1 + \beta_2 X_2 = 0\]

Any \(X = (X_1, X_2)^T\) for which this equation holds is a point on the hyperplane (line). This functional form generalizes to \(p\) dimensions quite easily:

\[\beta_0 + \beta_1 X_1 + \beta_2 X_2 + \dots + \beta_p X_p = 0\]

Again, for any point \(X = (X_1, X_2, \dots, X_p)^T\) in \(p\)-dimensional space (i.e. a vector of length \(p\)) that equals 0, then \(X\) lies on the hyperplane.

For \(X\) that does not meet this condition, then the data point lies on either side of the hyperplane:

\[\beta_0 + \beta_1 X_1 + \beta_2 X_2 + \dots + \beta_p X_p > 0\]

\[\beta_0 + \beta_1 X_1 + \beta_2 X_2 + \dots + \beta_p X_p < 0\]

The hyperplane therefore divides the \(p\)-dimensional space into two halves. To determine on which side of the hyperplane an observation lies, we simply calculate the sign of the corresponding hyperplane equation.

sim_hyper <- data_frame(x1 = seq(-1.5, 1.5, length.out = 20),
                        x2 = seq(-1.5, 1.5, length.out = 20)) %>%
  expand(x1, x2) %>%
  mutate(y = 1 + 2 * x1 + 3 * x2,
         group = ifelse(y < 0, -1,
                        ifelse(y > 0, 1, 0)),
         group = factor(group))

sim_hyper_line <- data_frame(x1 = seq(-1.5, 1.5, length.out = 20),
                             x2 = (-1 - 2 * x1) / 3)

ggplot(sim_hyper, aes(x1, x2, color = group)) +
  geom_point() +
  geom_line(data = sim_hyper_line, aes(color = NULL)) +
  labs(title = "Hyperplane in two dimensions") +
  theme(legend.position = "none")

Classification using a separating hyperplane

Let’s represent a hypothetical classification problem as the following: suppose we have an \(n \times p\) data matrix \(\mathbf{X}\) that consists of \(n\) training observations with \(p\) predictors in \(p\)-dimensional space:

\[x_1 = \begin{pmatrix} x_{11} \\ \vdots \\ x_{1p} \end{pmatrix}, \dots, x_n = \begin{pmatrix} x_{n1} \\ \vdots \\ x_{np} \end{pmatrix}\]

These observations fall into one of two classes \(y_1, \dots, y_n \in \{-1, 1 \}\) where \(-1\) and \(1\) represent two separate classes or categories. We also have a test observation \(x^*\) which is a \(p\)-vector of observed predictors \(x^* = (x_1^*, \dots, x_p^*)\). We want to develop a model that classifies the test observation correctly given our knowledge of the training observations. Previously we have used methods such as logistic regression (where the response variable is coded \(\{0, 1 \}\)) and decision trees to perform this task. Now we want to use a hyperplane to separate the training observations into the two possible classes.

A separating hyperplane perfectly separates training observations into their class labels. Observations in the blue class are coded as \(y_i = 1\) those from the red class as \(y_i = -1\). So a separating hyperplane takes on the properties:

\[\beta_0 + \beta_1 x_{i1} + \dots + \beta_p x_{ip} > 0, \text{if } y_i = 1\] \[\beta_0 + \beta_1 x_{i1} + \dots + \beta_p x_{ip} < 0, \text{if } y_i = -1\]

sim <- data_frame(x1 = runif(20, -2, 2),
                  x2 = runif(20, -2, 2),
                  y = ifelse(1 + 2 * x1 + 3 * x2 < 0, -1, 1)) %>%
  mutate_each(funs(ifelse(y == 1, . + 1.5, .)), x2) %>%
  mutate(y = factor(y, levels = c(-1, 1))) %>%
  mutate(line1 = (-1 - 2 * x1) / 3,
         line2 = .5 + (-1 - 1.5 * x1) / 3,
         line3 = .25 - .05 * x1)

ggplot(sim, aes(x1)) +
  geom_point(aes(y = x2, color = y)) +
  geom_line(aes(y = line1, color = NULL)) +
  geom_line(aes(y = line2, color = NULL)) +
  geom_line(aes(y = line3, color = NULL)) +
  labs(title = "Examples of separating hyperplanes") +
  theme(legend.position = "none")

If a separating hyperplane exists, then we can classify test observations based on their location relative to the hyperplane:

sim_mod <- svm(y ~ x1 + x2, data = sim, kernel = "linear", cost = 1e05,
               scale = FALSE)
sim_coef <- c(sim_mod$rho, t(sim_mod$coefs) %*% sim_mod$SV)

sim_grid <- data_frame(x1 = seq(-2, 2, length.out = 100),
                  x2 = seq(-2, 3.5, length.out = 100)) %>%
  expand(x1, x2) %>%
  mutate(y = ifelse(-sim_coef[[1]] + sim_coef[[2]] * x1 + sim_coef[[3]] * x2 > 0, -1, 1),
         y = factor(y, levels = c(-1, 1)))

sim_plane <- data_frame(x1 = seq(-2, 2, length.out = 100),
                        x2 = (sim_coef[[1]] - sim_coef[[2]] * x1) / sim_coef[[3]])

ggplot(sim, aes(x1)) +
  geom_point(data = sim_grid, aes(x1, x2, color = y), alpha = .25, size = .25) +
  geom_point(aes(y = x2, color = y)) +
  geom_line(data = sim_plane, aes(x1, x2)) +
  labs(title = "Maximal margin classification") +
  theme(legend.position = "none")

Classifications are based off the sign of \(f(x^*) = \beta_0 + \beta_1 x_1^* + \dots + \beta_p x_p^*\). If \(f(x^*)\) is positive, then we predict the test observation is \(1\). If \(f(x^*)\) is negative, then we predict the test observation is \(-1\). We can also consider the magnitude of \(f(x^*)\): the farther the magnitude is away from zero, then the farther the test observation falls from the hyperplane. We can be more confident of our predictions for observations far from the hyperplane, and less so for observations near the hyperplane (i.e. \(f(x^*)\) close to zero). The classifier resulting from the separating hyperplane \(f(x^*) = \beta_0 + \beta_1 x_1^* + \dots + \beta_p x_p^*\) is a linear decision boundary because the function itself is a linear form.

Maximal margin classifier

As we saw previously, if the data can be perfectly separated by a hyperplane it is likely true that there are multiple potential separating hyperplanes. We need a method for identifying the optimal separating hyperplane. This is known as the maximal margin hyperplane, which is the separating hyperplane that is farthest from the training observations. The margin is the smallest possible (perpendicular) distance between a training observation and the separating hyperplane. This distance is simply \(\hat{f}(x_i)\). The maximal margin hyperplane defines the hyperplane that minimizes the marginal distance across all training observations, and can be used to classify the test observation \(x^*\) based on which side of the hyperplane it lies. This is known as the maximal margin classifier. The expectation is that a classifier with a large margin for the training observations will also have a large margin for the test observations, leading to accurate classifications. As with the other methods we have discussed so far, this is an assumption and it is still possible to overfit the training data using the maximal margin classifier.

sim_pred <- predict(sim_mod, sim, decision.values = TRUE)
sim_dist <- attr(sim_pred, "decision.values")

ggplot(sim, aes(x1)) +
  geom_point(aes(y = x2, color = y)) +
  geom_point(data = sim_grid, aes(x1, x2, color = y), alpha = .1, size = .25) +
  geom_line(data = sim_plane, aes(x1, x2)) +
    geom_line(data = mutate(sim_plane, x2 = x2 - min(abs(sim_dist))),
              aes(x1, x2), linetype = 2) +
    geom_line(data = mutate(sim_plane, x2 = x2 + min(abs(sim_dist))),
              aes(x1, x2), linetype = 2) +
  labs(title = "Maximal margin classification") +
  theme(legend.position = "none")

Two observations are equidistant from the maximal margin hyperplane and lie along the dashed lines indicating the width of the margin. These observations are called the support vectors. They are vectors in \(p\)-dimensional space and “support” the maximal margin hyperplane because if the observations shifted at all in their predictor values \(X\), then the maximal margin hyperplane would shift as well. In fact, the maximal margin hyperplane is defined entirely by the support vectors; changes to the other observations would not effect the separating hyperplane as long as the changed observations do not cross the boundary set by the margin.

Constructing the maximal margin hyperplane

Constructing the maximal margin hyperplane is a (relatively) straight forward affair. Consider a set of \(n\) training observations with some number of real number predictors \(x_1, \dots, x_n \in \mathbb{R}^p\) and associated class labels \(y_1, \dots, y_n \in \{-1, 1\}\). We want to solve the optimization problem:

\[\begin{aligned} & \underset{\beta_0, \beta_1, \dots, \beta_p}{\text{maximize}} & & M \\ & \text{s.t.} & & \sum_{j=1}^p \beta_j^2 = 1, \\ & & & y_i(\beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + \dots + \beta_p x_{ip}) \geq M \; \forall \; i = 1, \dots, n \\ \end{aligned}\]

This is simpler than it looks. \(y_i(\beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + \dots + \beta_p x_{ip}) \geq M \; \forall \; i = 1, \dots, n\) requires the maximal margin hyperplane to sort observations on the correct side of the hyperplane with some amount of cushion, provided \(M\) is positive. The requirement \(\sum_{j=1}^p \beta_j^2 = 1\) means that not only are the observations sorted onto the correct sides of the hyperplane, but that the function \(y_i(\beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + \dots + \beta_p x_{ip})\) defines the perpendicular distance between the observation \(y_i\) and the hyperplane. Therefore \(M\) defines the margin of the hyperplane (i.e. the amount of cushion between the hyperplane and the closest training observations), so we select values for the parameters \(\beta_0, \beta_1, \dots, \beta_p\) to maximize \(M\); that is, obtain the largest amount of cushion possible given the training observations.

Non-separable cases

Unfortunately the maximal margin classifier only works if there exists a separating hyperplane for the data. If the cases cannot be perfectly separated by a hyperplane, then we can never satisfy the conditions of the maximal margin classifier.

data_frame(x1 = runif(20, -2, 2),
           x2 = runif(20, -2, 2),
           y = c(rep(-1, 10), rep(1, 10))) %>%
  mutate(y = factor(y, levels = c(-1, 1))) %>%
  ggplot(aes(x1, x2, color = y)) +
  geom_point() +
  labs(title = "Non-separable data") +
  theme(legend.position = "none")

Support vector classifier

Support vector classifiers relax the requirement of the maximal margin classifier by allowing the separating hyperplane to not perfectly separate the observations; instead, it can make some errors. This is reasonable when:

1. There exists no perfectly separating hyperplane
2. A perfectly separating hyperplane is too sensitive to individual training observations, generating potentially very small margins or overfitting the training set.²

# original model
sensitive <- data_frame(x1 = runif(20, -2, 2),
                  x2 = runif(20, -2, 2),
                  y = ifelse(1 + 2 * x1 + 3 * x2 < 0, -1, 1)) %>%
  mutate_each(funs(ifelse(y == 1, . + .5, .)), x2) %>%
  mutate(y = factor(y, levels = c(-1, 1)))

sens_mod <- svm(y ~ x1 + x2, data = sensitive, kernel = "linear",
                cost = 1e05, scale = FALSE)
sens_coef <- c(sens_mod$rho, t(sens_mod$coefs) %*% sens_mod$SV)
sens_plane <- data_frame(x1 = seq(-2, 2, length.out = 100),
                        x2 = (sens_coef[[1]] - sens_coef[[2]] * x1) / sens_coef[[3]])

ggplot(sensitive, aes(x1)) +
  geom_point(aes(y = x2, color = y)) +
  geom_line(data = sens_plane, aes(x1, x2)) +
  labs(title = "Maximal margin classification") +
  theme(legend.position = "none")

# slight tweak
sensitive2 <- data_frame(x1 = with(sensitive, x1[which(x2 == max(x2[y == -1]))]),
                         x2 = with(sensitive, max(x2[y == -1])) + .1,
                         y = factor(1, levels = c(-1, 1))) %>%
  bind_rows(sensitive)

sens2_mod <- svm(y ~ x1 + x2, data = sensitive2, kernel = "linear",
                cost = 1e05, scale = FALSE)
sens2_coef <- c(sens2_mod$rho, t(sens2_mod$coefs) %*% sens2_mod$SV)
sens2_plane <- data_frame(x1 = seq(-2, 2, length.out = 100),
                        x2 = (sens2_coef[[1]] - sens2_coef[[2]] * x1) / sens2_coef[[3]])

ggplot(sensitive2, aes(x1)) +
  geom_point(aes(y = x2, color = y)) +
  geom_line(data = sens2_plane, aes(x1, x2)) +
  geom_line(data = sens_plane, aes(x1, x2), linetype = 2) +
  labs(title = "Maximal margin classification") +
  theme(legend.position = "none")

Instead, we want a separating hyperplane that does not perfectly separate the two classes but provides greater robustness to individual observations and better classification of most training observations. We are willing to sacrifice accuracy on a few observations if the resulting hyperplane performs better across the remaining observations.

This approach is called the support vector classifier. It allows observations to not only exist on the wrong side of the margin (i.e. inside the cushion defined by \(M\)), but also exist on the wrong side of the hyperplane.

The approach is the same as the maximal margin classifier but the optimization problem is slightly different:

\[\begin{aligned} & \underset{\beta_0, \beta_1, \dots, \beta_p, \epsilon_1, \dots, \epsilon_n}{\text{maximize}} & & M \\ & \text{s.t.} & & \sum_{j=1}^p \beta_j^2 = 1, \\ & & & y_i(\beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + \dots + \beta_p x_{ip}) \geq M(1 - \epsilon_i), \\ & & & \epsilon_i \geq 0, \sum_{i = 1}^n \epsilon_i \leq C \\ \end{aligned}\]

As in the maximal margin classifier, we attempt to optimize \(M\) to generate the largest possible margin. However now we allow some error \(\epsilon_i\) for each observation so that they can fall on the wrong side of the margin or hyperplane.

• If \(\epsilon_i = 0\), then the \(i\)th observation falls on the correct side of the margin.
• If \(\epsilon_i > 0\), then the \(i\)th observation falls on the wrong side of the margin.
• If \(\epsilon_i > 1\), then the \(i\)th observation falls on the wrong side of the hyperplane.

\(C\) defines precisely how much error we are willing to tolerate in the resulting separating hyperplane. The sum of the errors for all training observations cannot exceed \(C\). Larger values of \(C\) permit more overall error in the separating hyperplane and lead to larger margins, and smaller values of \(C\) tolerate less error and produce smaller margins. If \(C = 0\) then we do not tolerate any error in the separating hyperplane, in which case \(\epsilon_1, \dots, \epsilon_n = 0\) and we estimate the maximal margin classifier (of course this is only possible if the classes are perfectly separable). Once we solve the optimization problem, we generate predictions the same way as for maximal margin classifiers, based on \(f(x^*) = \beta_0 + \beta_1 x_1^* + \dots + \beta_p x_p^*\).

Selecting a value for \(C\) is tricky and generally determined through a cross-validation approach to compare support vector classifiers under different values for \(C\). When \(C\) is small, we generate a model with low-bias (it fits the data well) but high-variance (small changes in the training observations can generate substantial changes in the support vector classifier). If \(C\) is large, we generate a model with more bias but less variance.

The important thing to realize is that the support vector classifier is robust, like the maximal margin classifier, to changes in observations outside of the margin. Observations that lie directly on the margin or inside the margin but on the correct side of the hyperplane are support vectors. The support vector classifier will only change if those observations are adjusted. When \(C\) is large, the number of observations falling inside the margin increases and therefore the number of support vectors also increases.

sim_c <- data_frame(x1 = rnorm(20),
                    x2 = rnorm(20),
                    y = ifelse(2 * x1 + x2 + rnorm(20, 0, .25) < 0, -1, 1)) %>%
  mutate(y = factor(y, levels = c(-1, 1)))

plot_svm <- function(df, cost = 1){
  # estimate model
  sim_mod <- svm(y ~ x1 + x2, data = df, kernel = "linear",
                 cost = cost,
                 scale = FALSE)
  
  # extract separating hyperplane
  sim_coef <- c(sim_mod$rho, t(sim_mod$coefs) %*% sim_mod$SV)
  sim_plane <- data_frame(x1 = seq(min(df$x1), max(df$x1), length.out = 100),
                          x2 = (-sim_coef[[1]] - sim_coef[[2]] * x1) / sim_coef[[3]])
  
  # extract properties to draw margins
  sim_pred <- predict(sim_mod, df, decision.values = TRUE)
  sim_dist <- attr(sim_pred, "decision.values")
  
  ggplot(df, aes(x1)) +
    geom_point(aes(y = x2, color = y)) +
    geom_line(data = sim_plane, aes(x1, x2)) +
    geom_line(data = mutate(sim_plane, x2 = x2 - min(abs(sim_dist))),
              aes(x1, x2), linetype = 2) +
    geom_line(data = mutate(sim_plane, x2 = x2 + min(abs(sim_dist))),
              aes(x1, x2), linetype = 2) +
    labs(subtitle = str_c("Cost = ", cost)) +
    coord_equal(xlim = range(df$x1),
                    ylim = range(df$x2)) +
    theme(legend.position = "none")
}

grid.arrange(grobs = list(plot_svm(sim_c, cost = 1),
                  plot_svm(sim_c, cost = 10),
                  plot_svm(sim_c, cost = 100),
                  plot_svm(sim_c, cost = 200)), ncol = 2)

Support vector machines

Non-linear decision boundaries

So far we have only demonstrated the support vector classifier with a linear decision boundary. But as with linear regression, we also know there are methods of extending the linear framework to account for non-linear relationships. Consider the following relationship:

set.seed(1)
x <- matrix(rnorm(200 * 2), ncol = 2)
x[1:100, ] <- x[1:100, ] + 2
x[101:150, ] <- x[101:150, ] - 2
y <- c(rep(1, 150), rep(2, 50))
sim_nonlm <- data.frame(x = x, y = as.factor(y)) %>%
  as_tibble %>%
  rename(x1 = x.1,
         x2 = x.2)

radial_p <- ggplot(sim_nonlm, aes(x1, x2, color = y)) +
  geom_point() +
  theme(legend.position = "none")
radial_p

A support vector classifier with a linear decision boundary would perform very poorly on this data.

We could go the route we discussed before and relax the linearity assumption by adding quadratic or cubic terms to address the non-linearity. For instance, adding a quadratic term would change the optimization problem to using \(2p\) features:

\[X_1, X_1^2, X_2, X_2^2, \dots, X_p, X_p^2\]

And therefore the optimization problem becomes:

\[\begin{aligned} & \underset{\beta_0, \beta_{11}, \beta_{12}, \dots, \beta_{p1}, \beta_{p2}, \epsilon_1, \dots, \epsilon_n}{\text{maximize}} & & M \\ & \text{s.t.} & & y_i \left( \beta_0 + \sum_{j = 1}^p \beta_{j1} x_{ij} + \sum_{j = 1}^p \beta_{j2} x_{ij}^2 \right) \geq M(1 - \epsilon_i), \\ & & & \epsilon_i \geq 0, \sum_{i = 1}^n \epsilon_i \leq C, \sum_{j = 1}^p \sum_{k = 1}^2 \beta_{jk}^2 = 1 \\ \end{aligned}\]

The problem with this approach is that as you add polynomial terms (or interactions or splines) you increase the feature space used to generate the decision boundary and the separating hyperplane (i.e. the total number of predictors increases). Maximizing this optimization problem is already computationally intensive: if you continue to increase the number of features, computing the support vector classifier becomes much more difficult and inefficient, and may even become impossible.

Support vector machines

The support vector machine is an extension of the support vector classifier that enlarges the feature space by using kernels. Kernels are a computationally efficient method for extending the feature space to accomodate a non-linear decision boundary.

Computing the support vector classifier involves the inner products of the observations, rather than the observations themselves.³ The inner product of two \(r\)-length vectors \(a\) and \(b\) is defined as \(\langle a,b \rangle = \sum_{i = 1}^r a_i b_i\).

(x <- 1:5)

## [1] 1 2 3 4 5

(y <- 1:5)

## [1] 1 2 3 4 5

x %*% y

##      [,1]
## [1,]   55

So the inner product of two observations is:

\[\langle x_i, x_{i'} \rangle = \sum_{j = 1}^p x_{ij} x_{i'j}\]

The linear support vector can be written as:

\[f(x) = \beta_0 + \sum_{i = 1}^n \alpha_i \langle x, x_i \rangle\]

where there are \(n\) parameters \(\alpha_i, i = 1, \dots, n\), one per training observation. To estimate the parameters \(\alpha_1, \dots, \alpha_n, \beta_0\), we just need to calculate the inner products between all pairs of training observations. However for observations which are not also support vectors, \(\alpha_i\) is actually zero. So in fact, we only need to calculate the inner products for support vectors \(\mathbb{S}\) which reduces the complexity of this task:

\[f(x) = \beta_0 + \sum_{i \in \mathbb{S}} \alpha_i \langle x, x_i \rangle\]

Kernels

Now rather than using the actual inner product,

\[\langle x_i, x_{i'} \rangle = \sum_{j = 1}^p x_{ij} x_{i'j}\]

instead we can use a generalization of the inner product following some functional form \(K\) which we will call a kernel:

\[K(x_i, x_{i'})\]

A kernel calculates the similarity of two observations. For example,

\[K(x_i, x_{i'}) = \sum_{j = 1}^p x_{ij} x_{i'j}\]

generates the support vector classifier, also known as the linear kernel. Alternatively, we could use a different kernel function such as:

\[K(x_i, x_{i'}) = (1 + \sum_{j = 1}^p x_{ij} x_{i'j})^d\]

This is called the polynomial kernel of degree \(d\) where \(d\) is some positive integer. This will generate a much more flexible decision boundary, similar to how using a spline in linear regression generates a flexible, non-linear functional form. To use this kernel in a support vector classifier, the functional form becomes:

\[f(x) = \beta_0 + \sum_{i \in \mathbb{S}} \alpha_i K(x,x_i)\]

sim_nonlm <- data_frame(x1 = runif(100, -2, 2),
                  x2 = runif(100, -2, 2),
                  y = ifelse(x1 + x1^2 + x1^3 - x2 < 0 +
                               rnorm(100, 0, 1), -1, 1)) %>%
  mutate(y = factor(y, levels = c(-1, 1)))

ggplot(sim_nonlm, aes(x1, x2, color = y)) +
  geom_point() +
  theme(legend.position = "none")

svm(y ~ x1 + x2, data = sim_nonlm, kernel = "polynomial", scale = FALSE, cost = 1) %>%
  plot(sim_nonlm, x2 ~ x1)

Another choice is the radial kernel:

\[K(x_i, x_{i'}) = \exp(- \gamma \sum_{j=1}^p (x_{ij} - x_{i'j})^2)\]

where \(\gamma\) is some positive constant. Radial kernels work by localizing predictions for test observations based on their Euclidian distance to nearby training observations.

sim_rad_mod <- svm(y ~ x1 + x2, data = sim_nonlm,
                     kernel = "radial", cost = 5, scale = FALSE)

radial_p

plot(sim_rad_mod, sim_nonlm, x2 ~ x1)

Kernels are better to use for support vector machines than other non-linear approachs because they do not enlarge the feature space. That is, you need to compute \(K(x_i, x_{i'})\) for all \(\binom{n}{2}\) distinct pairs \(i, i'\), but \(p\) itself remains the same. You do not need to explicitly enlarge the feature space to accomplish this task. The total number of features/predictors/independent variables in the model remains the same, so you can more easily compute the SVM.

Applying and interpreting SVMs

SVMs are generally used for prediction models. They generate predicted classes for test observations and we can assess confidence in the model and overall model fit using standard metric. However SVMs are not good for conducting inference, since there are no easy methods for interpreting the relative importance and influence of individual predictors on the separating hyperplane. Regression coefficients are generally easy to interpret, and even tree-based methods have visual and statistical interpretations (variable importance plots) of the individual predictors. Generally SVMs are interpreted by assessing overall model fit and error rates, using a combination of cross-validation methods and visuals such as ROC curves.

Titanic

Let’s try this method out on our trusty Titanic dataset, using age and gender to predict survival. First we’ll split our dataset into training and test sets.

titanic <- titanic_train %>%
  as_tibble %>%
  select(-Name, -Ticket, -Cabin, -PassengerId) %>%
  mutate_each(funs(as.factor(.)), Survived, Pclass, Embarked) %>%
  na.omit

titanic_split <- resample_partition(titanic, p = c("test" = .3, "train" = .7))

Our first attempt will use a linear kernel (i.e. support vector classifier) and we’ll use 10-fold cross-validation to determine the optimal cost parameter \(C\).

titanic_tune <- tune(svm, Survived ~ Age + Fare, data = as_tibble(titanic_split$train),
                     kernel = "linear",
                     range = list(cost = c(.001, .01, .1, 1, 5, 10, 100)))
summary(titanic_tune)

## 
## Parameter tuning of 'svm':
## 
## - sampling method: 10-fold cross validation 
## 
## - best parameters:
##  cost
##     5
## 
## - best performance: 0.364 
## 
## - Detailed performance results:
##    cost error dispersion
## 1 1e-03 0.400     0.0686
## 2 1e-02 0.390     0.0779
## 3 1e-01 0.376     0.0729
## 4 1e+00 0.366     0.0760
## 5 5e+00 0.364     0.0735
## 6 1e+01 0.364     0.0735
## 7 1e+02 0.364     0.0735

\(C = 1\) produces the lower CV error rate, so let’s use that model for estimating model fit using a ROC curve.

titanic_best <- titanic_tune$best.model
summary(titanic_best)

## 
## Call:
## best.tune(method = svm, train.x = Survived ~ Age + Fare, data = as_tibble(titanic_split$train), 
##     ranges = list(cost = c(0.001, 0.01, 0.1, 1, 5, 10, 100)), 
##     kernel = "linear")
## 
## 
## Parameters:
##    SVM-Type:  C-classification 
##  SVM-Kernel:  linear 
##        cost:  5 
##       gamma:  0.5 
## 
## Number of Support Vectors:  372
## 
##  ( 186 186 )
## 
## 
## Number of Classes:  2 
## 
## Levels: 
##  0 1

# get predictions for test set
fitted <- predict(titanic_best, as_tibble(titanic_split$test), decision.values = TRUE) %>%
  attributes

roc_line <- roc(as_tibble(titanic_split$test)$Survived, fitted$decision.values)
plot(roc_line)

auc(roc_line)

## Area under the curve: 0.759

How does this compare to a polynomial kernel SVM?

titanic_poly_tune <- tune(svm, Survived ~ Age + Fare, data = as_tibble(titanic_split$train),
                     kernel = "polynomial",
                     range = list(cost = c(.001, .01, .1, 1, 5, 10, 100)))
summary(titanic_poly_tune)

## 
## Parameter tuning of 'svm':
## 
## - sampling method: 10-fold cross validation 
## 
## - best parameters:
##   cost
##  0.001
## 
## - best performance: 0.394 
## 
## - Detailed performance results:
##    cost error dispersion
## 1 1e-03 0.394     0.0844
## 2 1e-02 0.398     0.0872
## 3 1e-01 0.396     0.0863
## 4 1e+00 0.398     0.0851
## 5 5e+00 0.396     0.0858
## 6 1e+01 0.398     0.0851
## 7 1e+02 0.396     0.0858

titanic_poly_best <- titanic_poly_tune$best.model
summary(titanic_poly_best)

## 
## Call:
## best.tune(method = svm, train.x = Survived ~ Age + Fare, data = as_tibble(titanic_split$train), 
##     ranges = list(cost = c(0.001, 0.01, 0.1, 1, 5, 10, 100)), 
##     kernel = "polynomial")
## 
## 
## Parameters:
##    SVM-Type:  C-classification 
##  SVM-Kernel:  polynomial 
##        cost:  0.001 
##      degree:  3 
##       gamma:  0.5 
##      coef.0:  0 
## 
## Number of Support Vectors:  397
## 
##  ( 198 199 )
## 
## 
## Number of Classes:  2 
## 
## Levels: 
##  0 1

# get predictions for test set
fitted <- predict(titanic_poly_best, as_tibble(titanic_split$test), decision.values = TRUE) %>%
  attributes

roc_poly <- roc(as_tibble(titanic_split$test)$Survived, fitted$decision.values)
plot(roc_poly)

auc(roc_poly)

## Area under the curve: 0.724

Not quite as good. The optimal cost parameter is smaller (\(.1\)), but the associated CV error rate is higher than the linear kernel and the resulting test AUC is smaller. How does this stack up against the radial kernel?

titanic_rad_tune <- tune(svm, Survived ~ Age + Fare, data = as_tibble(titanic_split$train),
                     kernel = "radial",
                     range = list(cost = c(.001, .01, .1, 1, 5, 10, 100)))
summary(titanic_rad_tune)

## 
## Parameter tuning of 'svm':
## 
## - sampling method: 10-fold cross validation 
## 
## - best parameters:
##  cost
##   100
## 
## - best performance: 0.326 
## 
## - Detailed performance results:
##    cost error dispersion
## 1 1e-03 0.400     0.0625
## 2 1e-02 0.400     0.0625
## 3 1e-01 0.358     0.0830
## 4 1e+00 0.342     0.0819
## 5 5e+00 0.330     0.0634
## 6 1e+01 0.328     0.0681
## 7 1e+02 0.326     0.0700

titanic_rad_best <- titanic_rad_tune$best.model
summary(titanic_rad_best)

## 
## Call:
## best.tune(method = svm, train.x = Survived ~ Age + Fare, data = as_tibble(titanic_split$train), 
##     ranges = list(cost = c(0.001, 0.01, 0.1, 1, 5, 10, 100)), 
##     kernel = "radial")
## 
## 
## Parameters:
##    SVM-Type:  C-classification 
##  SVM-Kernel:  radial 
##        cost:  100 
##       gamma:  0.5 
## 
## Number of Support Vectors:  334
## 
##  ( 165 169 )
## 
## 
## Number of Classes:  2 
## 
## Levels: 
##  0 1

# get predictions for test set
fitted <- predict(titanic_rad_best, as_tibble(titanic_split$test), decision.values = TRUE) %>%
  attributes

roc_rad <- roc(as_tibble(titanic_split$test)$Survived, fitted$decision.values)
plot(roc_rad)

auc(roc_rad)

## Area under the curve: 0.7

The radial improves upon both the polynomial and the linear SVMs. The CV error rate is lower and the test AUC is higher.

It’s easier to compare if we plot the ROC curves on the same plotting window:

plot(roc_line, print.auc = TRUE, col = "blue")
plot(roc_poly, print.auc = TRUE, col = "red", print.auc.y = .4, add = TRUE)
plot(roc_rad, print.auc = TRUE, col = "orange", print.auc.y = .3, add = TRUE)

Based on our predictions from the test set, the radial SVM performs the best on the AUC, followed by the linear SVM, and worst of all the polynomial SVM.

Voter turnout

Let’s test the SVM method on our voter turnout data. Again, let’s start by splitting the data into training and test sets.⁴

(mh <- read_csv("data/mental_health.csv") %>%
  mutate_each(funs(as.factor(.)), vote96, black, female, married) %>%
  na.omit)

## # A tibble: 1,165 × 8
##    vote96 mhealth_sum   age  educ  black female married inc10
##    <fctr>       <dbl> <dbl> <dbl> <fctr> <fctr>  <fctr> <dbl>
## 1       1           0    60    12      0      0       0  4.81
## 2       1           1    36    12      0      0       1  8.83
## 3       0           7    21    13      0      0       0  1.74
## 4       0           6    29    13      0      0       0 10.70
## 5       1           1    41    15      1      1       1  8.83
## 6       1           2    48    20      0      0       1  8.83
## 7       0           9    20    12      0      1       0  7.22
## 8       0          12    27    11      0      1       0  1.20
## 9       1           2    28    16      0      0       1  7.22
## 10      1           0    72    14      0      0       1  4.01
## # ... with 1,155 more rows

mh_split <- resample_partition(mh, p = c("test" = .3, "train" = .7))

SVM

Next let’s compare a few different SVM models. Again we’ll use 10-fold CV on the training set to determine the optimal cost parameter.

Linear kernel

mh_lin_tune <- tune(svm, vote96 ~ ., data = as_tibble(mh_split$train),
                    kernel = "linear",
                    range = list(cost = c(.001, .01, .1, 1, 5, 10, 100)))
summary(mh_lin_tune)

## 
## Parameter tuning of 'svm':
## 
## - sampling method: 10-fold cross validation 
## 
## - best parameters:
##  cost
##     5
## 
## - best performance: 0.294 
## 
## - Detailed performance results:
##    cost error dispersion
## 1 1e-03 0.315     0.0513
## 2 1e-02 0.315     0.0513
## 3 1e-01 0.300     0.0586
## 4 1e+00 0.296     0.0633
## 5 5e+00 0.294     0.0647
## 6 1e+01 0.294     0.0647
## 7 1e+02 0.294     0.0647

mh_lin <- mh_lin_tune$best.model
summary(mh_lin)

## 
## Call:
## best.tune(method = svm, train.x = vote96 ~ ., data = as_tibble(mh_split$train), 
##     ranges = list(cost = c(0.001, 0.01, 0.1, 1, 5, 10, 100)), 
##     kernel = "linear")
## 
## 
## Parameters:
##    SVM-Type:  C-classification 
##  SVM-Kernel:  linear 
##        cost:  5 
##       gamma:  0.125 
## 
## Number of Support Vectors:  513
## 
##  ( 258 255 )
## 
## 
## Number of Classes:  2 
## 
## Levels: 
##  0 1

fitted <- predict(mh_lin, as_tibble(mh_split$test), decision.values = TRUE) %>%
  attributes

roc_line <- roc(as_tibble(mh_split$test)$vote96, fitted$decision.values)
plot(roc_line)

auc(roc_line)

## Area under the curve: 0.777

Polynomial kernel

mh_poly_tune <- tune(svm, vote96 ~ ., data = as_tibble(mh_split$train),
                    kernel = "polynomial",
                    range = list(cost = c(.001, .01, .1, 1, 5, 10, 100)))
summary(mh_poly_tune)

## 
## Parameter tuning of 'svm':
## 
## - sampling method: 10-fold cross validation 
## 
## - best parameters:
##  cost
##     1
## 
## - best performance: 0.284 
## 
## - Detailed performance results:
##    cost error dispersion
## 1 1e-03 0.315     0.0389
## 2 1e-02 0.315     0.0389
## 3 1e-01 0.304     0.0388
## 4 1e+00 0.284     0.0459
## 5 5e+00 0.300     0.0397
## 6 1e+01 0.301     0.0309
## 7 1e+02 0.311     0.0339

mh_poly <- mh_poly_tune$best.model
summary(mh_poly)

## 
## Call:
## best.tune(method = svm, train.x = vote96 ~ ., data = as_tibble(mh_split$train), 
##     ranges = list(cost = c(0.001, 0.01, 0.1, 1, 5, 10, 100)), 
##     kernel = "polynomial")
## 
## 
## Parameters:
##    SVM-Type:  C-classification 
##  SVM-Kernel:  polynomial 
##        cost:  1 
##      degree:  3 
##       gamma:  0.125 
##      coef.0:  0 
## 
## Number of Support Vectors:  506
## 
##  ( 263 243 )
## 
## 
## Number of Classes:  2 
## 
## Levels: 
##  0 1

fitted <- predict(mh_poly, as_tibble(mh_split$test), decision.values = TRUE) %>%
  attributes

roc_poly <- roc(as_tibble(mh_split$test)$vote96, fitted$decision.values)
plot(roc_poly)

auc(roc_poly)

## Area under the curve: 0.73

Radial kernel

mh_rad_tune <- tune(svm, vote96 ~ ., data = as_tibble(mh_split$train),
                    kernel = "radial",
                    range = list(cost = c(.001, .01, .1, 1, 5, 10, 100)))
summary(mh_rad_tune)

## 
## Parameter tuning of 'svm':
## 
## - sampling method: 10-fold cross validation 
## 
## - best parameters:
##  cost
##     1
## 
## - best performance: 0.288 
## 
## - Detailed performance results:
##    cost error dispersion
## 1 1e-03 0.315     0.0535
## 2 1e-02 0.315     0.0535
## 3 1e-01 0.314     0.0509
## 4 1e+00 0.288     0.0524
## 5 5e+00 0.291     0.0561
## 6 1e+01 0.293     0.0568
## 7 1e+02 0.313     0.0537

mh_rad <- mh_rad_tune$best.model
summary(mh_rad)

## 
## Call:
## best.tune(method = svm, train.x = vote96 ~ ., data = as_tibble(mh_split$train), 
##     ranges = list(cost = c(0.001, 0.01, 0.1, 1, 5, 10, 100)), 
##     kernel = "radial")
## 
## 
## Parameters:
##    SVM-Type:  C-classification 
##  SVM-Kernel:  radial 
##        cost:  1 
##       gamma:  0.125 
## 
## Number of Support Vectors:  510
## 
##  ( 266 244 )
## 
## 
## Number of Classes:  2 
## 
## Levels: 
##  0 1

fitted <- predict(mh_rad, as_tibble(mh_split$test), decision.values = TRUE) %>%
  attributes

roc_rad <- roc(as_tibble(mh_split$test)$vote96, fitted$decision.values)
plot(roc_rad)

auc(roc_rad)

## Area under the curve: 0.749

plot(roc_line, print.auc = TRUE, col = "blue")
plot(roc_poly, print.auc = TRUE, col = "red", print.auc.y = .4, add = TRUE)
plot(roc_rad, print.auc = TRUE, col = "orange", print.auc.y = .3, add = TRUE)

SVM kernel	CV training error rate
Linear	0.294
Polynomial	0.284
Radial	0.288

This time the SVM with the highest AUC is the linear model, followed by the radial and then the polynomial SVM. Interestingly, the linear SVM had the highest training error rate (cross-validated), followed by radial, and then polynomial with the lowest error rate. These are cross-validated measures, so it’s not as if they should be heavily biased. However they are all within 1 percentage point of each other, so the differences may not actually be that substantial. Further exploration could be warranted here.

We could tinker with the parameters for the polynomial and radial kernel SVMs, adjusting the number of degrees in the polynomial SVM and testing different constants \(\gamma\) for the radial SVM, again using 10-fold CV to select the optimal values. Instead though, let’s see how the SVM with the highest AUC (linear) stacks up with some of the other statistical learning methods we could apply.

Logistic regression

mh_logit <- glm(vote96 ~ ., data = as_tibble(mh_split$train), family = binomial)
summary(mh_logit)

## 
## Call:
## glm(formula = vote96 ~ ., family = binomial, data = as_tibble(mh_split$train))
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -2.508  -1.049   0.534   0.855   1.948  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -4.00019    0.60162   -6.65  3.0e-11 ***
## mhealth_sum -0.08478    0.02795   -3.03   0.0024 ** 
## age          0.04332    0.00585    7.41  1.3e-13 ***
## educ         0.20305    0.03439    5.90  3.5e-09 ***
## black1       0.11542    0.23625    0.49   0.6252    
## female1      0.20224    0.16607    1.22   0.2233    
## married1     0.21052    0.18544    1.14   0.2563    
## inc10        0.06051    0.03143    1.93   0.0542 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1016.74  on 815  degrees of freedom
## Residual deviance:  876.12  on 808  degrees of freedom
## AIC: 892.1
## 
## Number of Fisher Scoring iterations: 4

fitted <- predict(mh_logit, as_tibble(mh_split$test), type = "response")
logit_err <- mean(as_tibble(mh_split$test)$vote96 != round(fitted))

roc_logit <- roc(as_tibble(mh_split$test)$vote96, fitted)
plot(roc_logit)

auc(roc_logit)

## Area under the curve: 0.778

The test error rate for the logistic regression model is 0.261.

Decision tree

mh_tree <- tree(vote96 ~ ., data = as_tibble(mh_split$train))
mh_tree

## node), split, n, deviance, yval, (yprob)
##       * denotes terminal node
## 
##  1) root 816 1000 1 ( 0 1 )  
##    2) age < 48.5 516  700 1 ( 0 1 )  
##      4) educ < 12.5 191  300 0 ( 1 0 )  
##        8) educ < 11.5 55   60 0 ( 1 0 ) *
##        9) educ > 11.5 136  200 1 ( 0 1 ) *
##      5) educ > 12.5 325  400 1 ( 0 1 )  
##       10) mhealth_sum < 4.5 259  300 1 ( 0 1 ) *
##       11) mhealth_sum > 4.5 66   90 0 ( 1 0 ) *
##    3) age > 48.5 300  300 1 ( 0 1 )  
##      6) educ < 12.5 153  200 1 ( 0 1 )  
##       12) inc10 < 1.08335 43   60 1 ( 0 1 ) *
##       13) inc10 > 1.08335 110  100 1 ( 0 1 ) *
##      7) educ > 12.5 147   60 1 ( 0 1 ) *

plot(mh_tree)
text(mh_tree, pretty = 0)

fitted <- predict(mh_tree, as_tibble(mh_split$test), type = "class")
tree_err <- mean(as_tibble(mh_split$test)$vote96 != fitted)

roc_tree <- roc(as.numeric(as_tibble(mh_split$test)$vote96), as.numeric(fitted))
plot(roc_tree)

auc(roc_tree)

## Area under the curve: 0.576

The test error rate for the decision tree model is 0.315.

Bagging

mh_bag <- randomForest(vote96 ~ ., data = as_tibble(mh_split$train),
                         mtry = 7)
mh_bag

## 
## Call:
##  randomForest(formula = vote96 ~ ., data = as_tibble(mh_split$train),      mtry = 7) 
##                Type of random forest: classification
##                      Number of trees: 500
## No. of variables tried at each split: 7
## 
##         OOB estimate of  error rate: 32.4%
## Confusion matrix:
##    0   1 class.error
## 0 89 168       0.654
## 1 96 463       0.172

varImpPlot(mh_bag)

fitted <- predict(mh_bag, as_tibble(mh_split$test), type = "prob")[,2]

roc_bag <- roc(as_tibble(mh_split$test)$vote96, fitted)
plot(roc_bag)

auc(roc_bag)

## Area under the curve: 0.714

Random forest

mh_rf <- randomForest(vote96 ~ ., data = as_tibble(mh_split$train))
mh_rf

## 
## Call:
##  randomForest(formula = vote96 ~ ., data = as_tibble(mh_split$train)) 
##                Type of random forest: classification
##                      Number of trees: 500
## No. of variables tried at each split: 2
## 
##         OOB estimate of  error rate: 30.4%
## Confusion matrix:
##    0   1 class.error
## 0 83 174       0.677
## 1 74 485       0.132

varImpPlot(mh_rf)

fitted <- predict(mh_rf, as_tibble(mh_split$test), type = "prob")[,2]

roc_rf <- roc(as_tibble(mh_split$test)$vote96, fitted)
plot(roc_rf)

auc(roc_rf)

## Area under the curve: 0.743

Compare the ROC curves

plot(roc_poly, print.auc = TRUE, col = "blue", print.auc.x = .2)
plot(roc_logit, print.auc = TRUE, col = "red", print.auc.x = .2, print.auc.y = .4, add = TRUE)
plot(roc_tree, print.auc = TRUE, col = "orange", print.auc.x = .2, print.auc.y = .3, add = TRUE)
plot(roc_bag, print.auc = TRUE, col = "green", print.auc.x = .2, print.auc.y = .2, add = TRUE)
plot(roc_rf, print.auc = TRUE, col = "purple", print.auc.x = .2, print.auc.y = .1, add = TRUE)

• SVM (linear kernel)
• Logistic regression
• Decision tree
• Bagging (\(n = 500\))
• Random forest (\(n = 500, m = \sqrt{p}\))

Based solely on the test AUC, logistic regression and random forest provides the highest predictive accuracy, slightly better than the linear kernel SVM. Decision tree performs the worst, though admittedly AUC is biased against it since all decision trees produce are predictions, not probabilities, so the ROC “curve” is actually a point.

Session Info

devtools::session_info()

##  setting  value                       
##  version  R version 3.5.1 (2018-07-02)
##  system   x86_64, darwin15.6.0        
##  ui       X11                         
##  language (EN)                        
##  collate  en_US.UTF-8                 
##  tz       America/Chicago             
##  date     2019-01-02                  
## 
##  package      * version date       source        
##  assertthat     0.2.0   2017-04-11 CRAN (R 3.5.0)
##  backports      1.1.2   2017-12-13 CRAN (R 3.5.0)
##  base         * 3.5.1   2018-07-05 local         
##  base64enc      0.1-3   2015-07-28 CRAN (R 3.5.0)
##  bindr          0.1.1   2018-03-13 CRAN (R 3.5.0)
##  bindrcpp       0.2.2   2018-03-29 CRAN (R 3.5.0)
##  broom        * 0.5.0   2018-07-17 CRAN (R 3.5.0)
##  cellranger     1.1.0   2016-07-27 CRAN (R 3.5.0)
##  class          7.3-14  2015-08-30 CRAN (R 3.5.1)
##  cli            1.0.0   2017-11-05 CRAN (R 3.5.0)
##  colorspace     1.3-2   2016-12-14 CRAN (R 3.5.0)
##  compiler       3.5.1   2018-07-05 local         
##  crayon         1.3.4   2017-09-16 CRAN (R 3.5.0)
##  datasets     * 3.5.1   2018-07-05 local         
##  devtools       1.13.6  2018-06-27 CRAN (R 3.5.0)
##  digest         0.6.18  2018-10-10 cran (@0.6.18)
##  dplyr        * 0.7.8   2018-11-10 cran (@0.7.8) 
##  e1071        * 1.7-0   2018-07-28 CRAN (R 3.5.0)
##  evaluate       0.11    2018-07-17 CRAN (R 3.5.0)
##  forcats      * 0.3.0   2018-02-19 CRAN (R 3.5.0)
##  gbm          * 2.1.3   2017-03-21 CRAN (R 3.5.0)
##  ggplot2      * 3.1.0   2018-10-25 cran (@3.1.0) 
##  glue           1.3.0   2018-07-17 CRAN (R 3.5.0)
##  graphics     * 3.5.1   2018-07-05 local         
##  grDevices    * 3.5.1   2018-07-05 local         
##  grid         * 3.5.1   2018-07-05 local         
##  gridExtra    * 2.3     2017-09-09 CRAN (R 3.5.0)
##  gtable         0.2.0   2016-02-26 CRAN (R 3.5.0)
##  haven          1.1.2   2018-06-27 CRAN (R 3.5.0)
##  hms            0.4.2   2018-03-10 CRAN (R 3.5.0)
##  htmltools      0.3.6   2017-04-28 CRAN (R 3.5.0)
##  httr           1.3.1   2017-08-20 CRAN (R 3.5.0)
##  ISLR         * 1.2     2017-10-20 CRAN (R 3.5.0)
##  jsonlite       1.5     2017-06-01 CRAN (R 3.5.0)
##  knitr          1.20    2018-02-20 CRAN (R 3.5.0)
##  lattice      * 0.20-35 2017-03-25 CRAN (R 3.5.1)
##  lazyeval       0.2.1   2017-10-29 CRAN (R 3.5.0)
##  lubridate      1.7.4   2018-04-11 CRAN (R 3.5.0)
##  magrittr       1.5     2014-11-22 CRAN (R 3.5.0)
##  Matrix         1.2-14  2018-04-13 CRAN (R 3.5.1)
##  memoise        1.1.0   2017-04-21 CRAN (R 3.5.0)
##  methods      * 3.5.1   2018-07-05 local         
##  modelr       * 0.1.2   2018-05-11 CRAN (R 3.5.0)
##  munsell        0.5.0   2018-06-12 CRAN (R 3.5.0)
##  nlme           3.1-137 2018-04-07 CRAN (R 3.5.1)
##  parallel     * 3.5.1   2018-07-05 local         
##  pillar         1.3.0   2018-07-14 CRAN (R 3.5.0)
##  pkgconfig      2.0.2   2018-08-16 CRAN (R 3.5.1)
##  plyr           1.8.4   2016-06-08 CRAN (R 3.5.0)
##  pROC         * 1.12.1  2018-05-06 CRAN (R 3.5.0)
##  purrr        * 0.2.5   2018-05-29 CRAN (R 3.5.0)
##  R6             2.3.0   2018-10-04 cran (@2.3.0) 
##  randomForest * 4.6-14  2018-03-25 CRAN (R 3.5.0)
##  rcfss        * 0.1.5   2018-05-30 local         
##  Rcpp           1.0.0   2018-11-07 cran (@1.0.0) 
##  readr        * 1.1.1   2017-05-16 CRAN (R 3.5.0)
##  readxl         1.1.0   2018-04-20 CRAN (R 3.5.0)
##  rlang          0.3.0.1 2018-10-25 CRAN (R 3.5.0)
##  rmarkdown      1.10    2018-06-11 CRAN (R 3.5.0)
##  rprojroot      1.3-2   2018-01-03 CRAN (R 3.5.0)
##  rstudioapi     0.7     2017-09-07 CRAN (R 3.5.0)
##  rvest          0.3.2   2016-06-17 CRAN (R 3.5.0)
##  scales         1.0.0   2018-08-09 CRAN (R 3.5.0)
##  splines      * 3.5.1   2018-07-05 local         
##  stats        * 3.5.1   2018-07-05 local         
##  stringi        1.2.4   2018-07-20 CRAN (R 3.5.0)
##  stringr      * 1.3.1   2018-05-10 CRAN (R 3.5.0)
##  survival     * 2.42-6  2018-07-13 CRAN (R 3.5.0)
##  tibble       * 1.4.2   2018-01-22 CRAN (R 3.5.0)
##  tidyr        * 0.8.1   2018-05-18 CRAN (R 3.5.0)
##  tidyselect     0.2.5   2018-10-11 cran (@0.2.5) 
##  tidyverse    * 1.2.1   2017-11-14 CRAN (R 3.5.0)
##  titanic      * 0.1.0   2015-08-31 CRAN (R 3.5.0)
##  tools          3.5.1   2018-07-05 local         
##  tree         * 1.0-39  2018-03-17 CRAN (R 3.5.0)
##  utils        * 3.5.1   2018-07-05 local         
##  withr          2.1.2   2018-03-15 CRAN (R 3.5.0)
##  xml2           1.2.0   2018-01-24 CRAN (R 3.5.0)
##  yaml           2.2.0   2018-07-25 CRAN (R 3.5.0)

---

1. Though they can also be applied to regression on continuous response variables.↩

2. Remember that we can use the perpendicular distance from the hyperplane as a measure of confidence in our predictions, so the new training observation diminishes our confidence for quite a few of the red training observations.↩

3. Like how boosting uses the residuals of the response variable \(Y\), rather than \(Y\) itself.↩

4. Why use the validation set approach? We’ve discussed the inadequacies of it before. We could use \(k\)-fold cross validation instead, however setting this up with the proper code would be much more complicated. When conducting exploratory analysis, you don’t necessarily need to do this on your first pass through the data. Certainly I recommend using CV to validate your models and compare them before publishing anything, but for this application I think the validation set approach works fine.↩

This work is licensed under the CC BY-NC 4.0 Creative Commons License.